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Concepts of information theory are increasingly used to characterize collective phenomena in condensed 
matter systems, such as the use of entanglement entropies to identify emergent topological order in interact¬ 
ing quantum many-body systems. Here we employ classical variants of these concepts, in particular Renyi 
entropies and their associated mutual information, to identify topological order in classical systems. Like for 
their quantum counterparts, the presence of topological order can be identified in such classical systems via a 
universal, subleading contribution to the prevalent volume and boundary laws of the classical Renyi entropies. 
We demonstrate that an additional subleading 0(1) contribution genetically arises for all Renyi entropies 5^"' 
with n > 2 when driving the system towards a phase transition, e.g. into a conventionally ordered phase. This 
additional subleading term, which we dub connectivity contribution, tracks back to partial subsystem ordering 
and is proportional to the number of connected parts in a given bipartition. Notably, the Levin-Wen summation 
scheme - typically used to extract the topological contribution to the Renyi entropies - does not fully elim¬ 
inate this additional connectivity contribution in this classical context. This indicates that the distillation of 
topological order from Renyi entropies requires an additional level of scrutiny to distinguish topological from 
non-topological 0(1) contributions. This is also the case for quantum systems, for which we discuss which 
entropies are sensitive to these connectivity contributions. We showcase these findings by extensive numerical 
simulations of a classical variant of the toric code model, for which we study the stability of topological order 
in the presence of a magnetic field and at finite temperatures from a Renyi entropy perspective. 


I. INTRODUCTION 

The ground states of interacting many-body systems can ex¬ 
hibit subtle forms of order such as the formation of long-range 
topological order or the emergence of spin liquid physics^ - 
both in the quantum and classical world. Distilling the precise 
nature of these unconventional forms of order or even identi¬ 
fying their mere existence from, for instance, the ground state 
wavefunction of a quantum many-body system is a highly 
non-trivial endeavor. This is particularly true if one tries to 
rely on conventional characterization approaches (often in¬ 
spired by their experimental feasibility) such as the use of 
correlation functions. As it is oftentimes the case, the break¬ 
through to overcome these limitations has come by completely 
changing one’s point of view - in this case by adopting an in¬ 
formation theory perspective on the many-body systerrP. 

One core concept of quantum information theory is to con¬ 
sider the entanglement for a bipartition of the many-body sys¬ 
tem into two parts and to precisely map out its dependencies 
on the size and shape of the boundary separating the two parts. 
This idea has been pioneered by Bekenstein and Hawking 
who in their description of black holes have introduced the 
entanglement entropy as a quantitative measure of entangle¬ 
ment and demonstrated its characteristic boundary law scal¬ 
ing. The best-known example of an entanglement entropy is 
the von-Neumann entropy 

S'(A) =-Tr(p^lnp^) , (1) 

where pA is the reduced density matrix of subsystem A after 
tracing out subsystem B. By considering higher moments of 
the density matrix, one can embed the von-Neumann entropy 
into a family of Renyi entropies 

Sn(.A) = ^\n{Tr[p'X]) , ( 2 ) 


where the index n typically is an integer n> 2 with the limit 
n —^ 1 recovering the von-Neumann entropy Q. 

In contrast to thermal entropies these entanglement en¬ 
tropies are not extensive, but show a prevalent scaling with 
the length i of the boundary of the bipartitiorP 

Sn{A) = a ■ i + b-Ini — ^ + 0{l/i), (3) 

where in addition to the boundary law scaling of the first term 
we explicitly list the first two subleading terms. The loga¬ 
rithmic corrections typically arise from corner contributions 
for gapless phases or quantum critical pointP. We will fo¬ 
cus in the following on the 0(1) contribution 7 that is typi- 
cally associated with the topological entanglement entropjEMl 
for gapped two-dimensional phases, though such 0(1) contri¬ 
butions are also known to arise at quantum critical points re¬ 
flecting geometrical aspects of the subsysterrP. For a gapped 
system a finite value of 7 not only indicates the presence of 
topological order, but also allows to narrow dowrP the pre¬ 
cise type of the underlying topological quantum held theory 
due to its universal character^. One way to calculate the en¬ 
tanglement entropy 7 from Equation 0 is to perform a care¬ 
ful scaling analysis for bipartitions of varying size, as it was 
done, for instance, in identifying the emergence of topological 
order in the ground state of the Heisenberg antiferromagnet 
on the kagome lattic^S^. Alternatively, one can compute the 
topological entanglement entropy using one of the summa¬ 
tion schemes put forward by Levin and WerP or Kitaev and 
PreskilPthat consider a set of bipartitions of a system of fixed 
system size constructed in a way that all leading boundary and 
subleading corner contributions to the entanglement entropy 
precisely cancel. An illustration of the four bipartitions used 
in the Levin-Wen scheme is provided in Lig. 

A similar information theory perspective can also be taken 
to identify subtle orders in classical many-body systems, 
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FIG. 1. (color online) Illustration of the four bipartitions used in 
the Levin-Wen summation scheme to extract the topological entropy. 
The latter is calculated as 7 = —^(Ai)-|- 5 (^ 2 ) + •S'(^ 3 ) — 5 '(^ 4 ). 


which will be the subject of this manuscript. In making this 
quantum to classical transition, we follow the work of Castel- 
novo and ChamoiP^ and consider the Shannon entropy 

= - XI PoA InPjA : (4) 

{ja} 

where {Ja} denotes the set of possible configurations of sub¬ 
system A in a bipartition of the system with statistical weights 
Pj^. More generally, we can consider an entire family of clas¬ 
sical Renyi entropies 

\{U} / 

where n again is a positive integer n > 2 and the limes n —1 
recovers the Shannon entropy 0. 

One key distinction of the classical versus quantum Renyi 
entropies is that the classical Renyi entropies follow a volume 
law 


Sn{A) = cv ■Va + ce ■ £-J + 0{l/i), (6) 

where Va now indicates the volume of partition A. The sub¬ 
leading 0(1) contribution 7 again indicates the existence of 
topological order. For the classical toric code, which we will 
discuss in more detail in the following, its value of 7 = In 2 
resembles the result of the quantum version up to a factor of 2 
as first discussed in Ref.[TO]and recently expanded to a family 
of more general classical stringnet model^lM 

In this manuscript, we carefully analyze this 0(1) contri¬ 
bution to the classical Renyi entropies when driving a sys¬ 
tem through a phase transition. As one particular example 
we consider the continuous phase transition from a topolog¬ 
ically ordered to a conventionally ordered phase, which in 
the classical toric code can be driven by a magnetic field h. 
One of our main results is that for all higher-order Renyi en¬ 
tropies with n > 2 there is an additional constant term 
beyond the topological contribution for the intermediate cou¬ 
pling regime hc/n < h < he- This additional contribution 
arises from partial subsystem ordering in this coupling regime. 
Because it is sensitive to the number of disconnected parts of 
these subsystems we dub it a “connectivity contribution” to 
sharply distinguish it from the known topological contribu¬ 
tion. To illustrate the generic, non-topological nature of this 
0(1) connectivity contribution we also demonstrate its occur¬ 
rence for the two-dimensional Ising model driven through its 
thermal phase transition. Finally, we come back to the classi¬ 
cal toric code and round off our discussion by analyzing the 


finite-temperature behavior of the 0(1) contributions to clas¬ 
sical Renyi entropies in a variety of settings. 

The manuscript is organized as follows. We start with an 
introduction to the classical toric code in Sec. [Ill We then 
turn to the information theory perspective of classical many- 
body systems and (re)introduce classical Renyi entropies and 
related measures such as the mutual information and discuss 
their scaling behavior in Sec. Ill In this section, we also mo¬ 
tivate and discuss our main results from an analytical perspec¬ 
tive. This is followed by a discussion of extensive numeri¬ 
cal simulations in Sec. EYl We conclude with an outlook on 
the applicability of our results to quantum many-body sys¬ 
tems and a discussion of the general shortcomings of Renyi 
entropies in Sec.|V] The main body of the manuscript is com¬ 
plemented by an extensive appendix, which gives a detailed 
exposition of the implementation of the replica technique to 
calculate the Renyi entropies. In particular, we comprehen¬ 
sively discuss how to combine non-local loop update tech¬ 
niques with the replica scheme. 


II. CLASSICAL TORIC CODE MODEL 


To set the stage for our analysis we give a short introduc¬ 
tion to the classical toric code model, which will serve as a 
paradigmatic example of a classical system with topological 
order in our subsequent discussion. The classical toric code is 
derived from its well known quantum counterpari^^, which is 
defined by the Hamiltonian 

4 4 

p^P i—1 v^V i—1 

where SU(2) spin-1/2 degrees of freedom located on the bonds 
of a square lattice are interacting via four-spin interactions 
around the plaquettes P and vertices V of the lattice. The 
two (commuting) terms energetically favor closed loop con¬ 
figurations either in the ax- or tr^-basis (which are dual to 
each other). The ground state of the Hamiltonian can then be 
written as a loop gas, i.e. an equal-weight superposition of all 
closed loop configurations (in one of the two bases). We will 
work in the a^-basis in the following. 

The classical toric code is defined by interpreting the loop 
gas as a classical partition function, i.e. a partition sum to 
which all closed loop configurations contribute equally. We 
can define a Hamiltonian for this classical system in terms of 
Ising variables located on the bonds of a square lattice as 

4 

H = ( 8 ) 

vGV i—1 

where a four-spin term for each vertex again favors even par¬ 
ity spin alignment corresponding to closed loop configura¬ 
tions. The low-temperature manifold of states for this clas¬ 
sical Hamiltonian is thus again providing us with a loop gas. 
Like in the quantum model we can define four distinct wind¬ 
ing number parity sectors, which for a toroidal geometry (pe¬ 
riodic boundary conditions) give rise to a 4-fold degeneracy at 
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zero temperature - a hallmark of topological order also in the 
classical context. In the following, we will initially focus on 
this zero-temperature physics and impose a hard-constraint on 
the system enforcing closed loop configurations. 

We can drive the classical toric code out of its topologically 
ordered phase by adding an external field term —h 
to the Hamiltonian that breaks the equal-weightness of the 
closed loop configurations constituting the loop gas. States 
with a higher total magnetization m = preferred 

via statistical weights (in lieu of uniform weights 1 for 
the unperturbed loop gas). Since a majority of spins pointing 
up in the toric code lattice effectively shortens the length of 
the individual loops, we will also call h a loop tension. For 
sufficiently large strength this loop tension will drive a tran¬ 
sition from the loop gas to a conventionally polarized state 
where all spins are aligned along the external field. 

Note that in defining these statistical weights we have not 
introduced the notion of temperature. We will return to this 
point in Sec. IV D where we discuss the role of open loop 
defects arising at finite temperatures in Hamiltonian ([^ and 
their interplay in the presence of a magnetic field. 


A. Kramers-Wannier duality and the 2D Ising model 

The zero-temperature physics of the classical toric code 
in a magnetic field can be mapped via a Kramers-Wannier 
duality^] to a 2D Ising model at finite temperatures - much 
like its quantum counterpart in a magnetic field can be mapped 
to a 3D Ising model at finite temperature^. To illustrate this 
mapping consider an alternative set of Ising degrees of free¬ 
dom Pi located at the center of the plaquettes P of the original 
toric code model, see Fig. For a given loop configuration 
{a} the values of these alternate Ising spins are chosen such 
that pipj = —where site i and j are nearest neighbors 
and cr(i,j) is the spin on the toric code bond that separates 
plaquettes i and j. The magnetic field term of the toric code 
model thus induces an Ising interaction between the dual de¬ 
grees of freedom. Note that loops in the toric code correspond 
to domain walls in the Ising model as illustrated in Fig. 

This duality between the toric code and the Ising model 
readily implies that the critical value of the loop tension he 
is directly related to the analytically knowrP^ critical value /3c 
of the finite-temperature transition in the Ising model 

/jTC ^ Rising ^ + V2)/2^ 0.440686.... 

It is also obvious that the mapping is unique up to a global 
flip of the Ising spins. In other words, the Z 2 symmetry of the 
Ising model is lost in the toric code model.G^ 

The mapping also lets us understand the phase transition in 
the classical toric code as a continuous phase transition (in the 
2D Ising universality class). In principle, one can of course 
track this transition by measuring, for instance, the total mag¬ 
netization m which serves as order parameter in the perspec¬ 
tive of the Ising model. However, here we want to pursue a 
different direction and ask how other global quantities such 
as the Renyi entropies can be used to quantitatively track this 
phase transition. 


Toric Code 


Ising model 



FIG. 2. (color online) The well-known Kramers-Wannier dualitjli^ 
maps the classical toric code in a magnetic field to the 2D Ising model 
at finite temperatures. 


III. CLASSICAL RENYI ENTROPIES 

With the classical toric code at hand, we now direct our 
analysis towards the signatures of topological order in clas¬ 
sical Renyi entropies. In this section our emphasis will be 
mostly on an analytical perspective. To this end, we first pro¬ 
vide a short introduction to the replica description of Renyi 
entropies. We then discuss the 0(1) contributions arising for 
the toric code model, in particular when driving the system 
through the field-induced phase transition to a topologically 
trivial phase. 


A. Replica technique 

Let us start our more detailed discussion of Renyi entropies 
by briefly laying out how one can calculate the Renyi en¬ 
tropies using the so-called replica techniqu^. The essen¬ 
tial idea of this approach is to calculate the n-th Renyi en¬ 
tropy from n replicas of the representation of the underlying 
many-body system. While bipartition B is allowed to inde¬ 
pendently fluctuate across the n replicas, the fluctuations of 
bipartition A are constrained to meet certain boundary condi¬ 
tions between the replicas. This scheme was original devel¬ 
oped in the context of analytical calculation^ of Renyi en¬ 
tropies for quantum held theories (on an n-sheeted Riemann 
surface) and later adapted to the analysis quantum many-body 
systems in condensed matter physic^. It is also directly 
amenable to numerical simulations and has first been ported 
to Monte Carlo simulations in the context of lattice gauge 
theorieJSl. For quantum lattice models the replica technique 
has first been reformulated in the context of stochastic series 
expansions in Ref. Subsequently the replica technique 
has been adapted to various other flavors of quantum Monte 
Carlo including variational Monte Carlcl^, continuum-space 
path integral Monte CarlcP^, determinantal Monte Carld^, 
fermionic continuous-time path integral Monte CarlcP^, and 
hybrid Monte CarlcP^. More recently implementations of 
this replica technique have been envisioned for experimental 
setting^, which has led to the first experimental measure¬ 
ment of Renyi entropies in an ensemble of ultracold atom^. 

Here we follow the worlsP^by laconis and collaborators and 
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adapt the replica technique to classical systems - an approach, 
which is again directly amenable to Monte Carlo simulations. 
We focus on the second Renyi entropy iS'2(4l) and therefore 
consider two replicas of the system where the degrees of free¬ 
dom in part A are precisely mimicking each other in the two 
replicas, while those in part B are allowed to fluctuate inde¬ 
pendently for the two replicas - see Fig. for an illustration. 
This replica representation of the Renyi entropy follows di¬ 
rectly from writing the second order Renyi entropy via its 
definition (|^ and employing statistical weights exp(/im) to 
capture the effect of a finite magnetic held of strength h 


52(41) = -In 


= - In 


= - In 


'y ^ (PctaUctb ) (PctaUctb/) 


^ A')^ B B' 

E 

Z[A,2,h] 

z[hV 


^h{TnA+mB+mA+'mgi) 


(9) 


Here we have introduced the partition function Z[h] of the 
total system in the presence of a magnetic held of strength h. 
The second partition function Z\A^ 2, h] captures a setting in 
which subsystem A is subject to a magnetic held of strength 
2 h, i.e. twice the strength of the applied magnetic held, while 
the two replicas for subsystem B are subject to a magnetic 
field of strength h. 

This replica representation of arbitrary Renyi entropies with 
n > 2 is then obtained in a straight-forward manner. In gen¬ 
eral, one encounters partition functions of the form Z[A, n, h] 
describing a system, in which subsystem A experiences a 
magnetic held of strength nh, while being coupled to n in¬ 
dependent replicas of subsystem B, each experiencing a mag¬ 
netic held of strength h. 



FIG. 3. Illustration of the replica representation of the second Renyi 
entropy 52(A). Pictured are two replicas of the toric code model 
with a possible loop gas configuration superimposed. Similar to 
Fig. 2 the thick black lines represent up-pointing spins on the bonds 
of the lattice, while the down-pointing spins are not illustrated. Note 
that the loop configuration in the connected part A is identical in both 
replicas, while two independent loop configurations are allowed to 
occur in the two replicas of part B. 


B. Topological entropy and connectivity contribution 

As mentioned already in the introduction, the classical 
Renyi entropies generic ally follow a volume law as given in 
Eq. @ augmented by a subleading boundary law and 0(1) 
contributions indicative of topological order. We will concen¬ 
trate on precisely these 0(1) contributions in the following. 
To extract them from the Renyi entropy we employ the Levin- 
Wen summation schemd^considering the four bipartitions Ai, 
A 2 , A 3 and A 4 illustrated in Fig.[2 Adding up the Renyi en¬ 
tropies of these four bipartitions as 

A5 = -52(41i) + 52 ( 412 ) + 52 (^ 3 ) - 52 (^ 4 ), (10) 

we note that all volume, boundary, and possible corner con¬ 
tributions to the Renyi entropies precisely cancel. Hence, A5 
provides us with a direct measurement of the remaining 0(1) 
contributions. Following the work of Ref. m we first note 
that in the classical context only the fourth bipartition in the 
Levin-Wen scheme contributes a non-zero topological contri¬ 
bution as listed in the Table of Fig. It is this non-vanishing 
contribution, which provides us with a topological entropy of 
precisely 7 = In 2 for the unperturbed classical toric code as 
originally shown by Castelnovo and ChamorPSI. 

We now consider the effect of a finite loop tension (or 
magnetic field). From the Kramers-Wannier duality of the 
Ising model we already know that only when the loop tension 
reaches the critical value he = ln(l -F s/ 2 ) 12 topological or¬ 
der in the system is destroyed. When considering the regime 
h < he from a renormalization group perspective, the effect 
of a loop tension renormalizes to zero in the thermodynamic 
limit and this physics should also be reflected in the entropy. 
One might thus be (mis)led to expect a constant 7 = In 2 con¬ 
tribution to the classical entropy all the way up to the critical 
value he- 

In fact, however, the situation turns out to be a bit more 
subtle for Renyi entropies. To get the topological entropy, it 
is mandatory that all replicas remain topologically ordered. 
But as we have seen in the replica representation of the Renyi 
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FIG. 4. (color online) Constant contributions to the second Renyi 
entropy in the topologically ordered phase of the classical toric code 
model in a magnetic field of strength h for the different bipartitions 
in the Levin-Wen scheme. For simplicity n = 2 is shown. 
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entropies in Sec. Ill A the two subsystems are exposed to dif¬ 
ferent strengths of the loop tension, e.g. for the second Renyi 
entropy one has field strengths 2 h and h for subsystems A 
and B, respectively. This requires to refine our (renormaliza¬ 
tion group) perspective to consider three separate parameter 
regions. 


When h < hc/^., all replicas are in the topological phase. 
Here the loop tension is an itTelevant perturbation (in an renor¬ 
malization group sense), and the arguments of Ref [TO] apply. 
Hence we expect a purely topological 0(1) contribution of 
7 = In 2 to be calculated by the Levin-Wen scheme. 


When h > he all replicas are in the ordered phase, with 
most spins pointing up. As can be understood by considering 
the limiting case h ^ oo with zero entropy, the four contri¬ 
bution of Eq. (lOi become identical and the 0(1) contribution 
AS” vanishes in the thermodynamic limit. 


The intermediate region he/2 < h < he is. somewhat more 
delicate. Here we have one replica A already in the ordered 
phase, while coupled to n replicas of B still in the topological 
phase. In case the bipartition leads to a single boundary (i.e. 
bipartitions A 2 and A 3 ) this generates an extra contribution 
In 2/(n — 1) to the entropy. In case the bipartition gives two 
separate boundaries (i.e. bipartitions Ai and A4, in which ei¬ 
ther subsystem A 01 B comes in two disconnected parts) the 
contribution becomes 2 In 2/(n — 1). We name this contribu¬ 
tion connectivity contribution, as it has no topological origin 
and is instead sensitive to the disconnected parts in subsys¬ 
tem Aot B. The various contributions to the entropy for each 
geometry are detailed in Table. specifying the situation of 
n = 2 for simplicity. 


Adding up the individual contributions, the Levin-Wen 
summation scheme then yields 


{ In 2 , h < he/2 

;^ln2 , he/2<h<he (11) 

0 , h>he 

in the three different parameter regions. 

This is one of the main results of this paper. It will be 
checked numerically for n = 2 and large system sizes in 
Sec. HYl Let us finally comment on the two special values 
h = he /2 and h — he- In that case, one or several of the 
replicas are exactly at the critical temperature in the language 
of the dual Ising model. These replicas are now governed by 
the Ising conformal field theory. It is knowrP^that at the crit¬ 
ical point the corner contributions are enhanced to logarith¬ 
mic terms cx ln(L/a), where a is some microscopic cutoff 
of the order of a lattice spacing. As is explained in Ref. 
these logarithmic terms are canceled by the linear combina¬ 
tion However, there is no reason for all microscopic cut¬ 
offs to be the same. This means there is going to be an extra 
non-universal 0(1) contribution to A5' arising in the vicin¬ 
ity of the critical points. Hence we expect a non-vanishing, 
but non-universal value for AS” exactly at the critical points 
h = he/2 and h = he¬ 


ll 11 II II I. I, I I. ( II I, I, 



FIG. 5. (color online) Fourth bipartition of the Levin Wen scheme as 
an example of how the subdivision is performed. Spins in subsystem 
A are shown in blue, spins in B are in red. The loops are represented 
by thick black lines. 


IV. NUMERICAL SIMULATIONS 


To complement the analytical perspective on Renyi en¬ 
tropies for the loop tension driven phase transition in the toric 
code we will now turn to a discussion of extensive numerical 
simulations of this model. 


A. Monte Carlo setup 


We start our discussion of the numerical simulations by first 
outlining some aspects of our classical Monte Carlo setup, in 
particular with regard to the choice of bipartition, the numer¬ 
ical adaptation of the replica technique as well as the imple¬ 
mentation of non-local loop update techniques. 


Bipartition 


Let us first discuss our general choice for bipartitioning the 
system. A naive approach would be to cut along the lines of 
the square lattice and therefore also cut through spins. This is 
certainly not desirable, as it would also imply that some loop 
segments could mn along the boundary of the bipartition (see 
also Fig.|^ and in general it would leave us wondering how to 
assign the degrees of freedom on the boundary to the respec¬ 
tive subsystems. Instead we choose to rotate the toric code 
lattice by 45 degrees and perform diagonal cuts in a way such 
that each plaquette is cut symmetrically along each line of the 
bipartition. This is illustrated for subsystem A 4 in Fig. 




6 


Replica technique 


We employ two alternative approaches to implement the 
replica representation (|^ and calculate the Renyi entropy 
from the ratio of the partition functions Z\h] and Z[A,2, h\ 
as 


S2{A) = - In 


Z[A,2,h] 

Z[h? 


in a classical Monte Carlo setup. 

The first approach uses the fact that for a given partition 
function, e.g. Z[h], we have = {m), where (m) is 

the average total magnetization. We can thus obtain the en¬ 
tropy by performing an integration over the estimated values 
of (m) for a dense range of parameters h. This method is an 
adaptation of the thermodynamic integ ration employed e.g. in 
the context of thermal phase transition^^lIII] 

The second approach is often referred to as ensemble 
switchin^^ as it relates the ratio of the two partition func¬ 
tions Z[hY and Z[A^ 2, h] to the number of times one can 
switch between the two ensembles for a given sequence of 
configurations in a Markov chain. To this end, one samples 
configurations that belong to the system with partition func¬ 
tion Z\hY, i.e. two copies of the full system. Since the con¬ 
figuration space of Z\h\’^ fully includes the one of Z\A,2, h\ 
we can readily estimate the relative weight of the two parti¬ 
tion functions (needed for the above ratio) by simply counting 
how often a configuration in the Markov chain of the parti¬ 
tion function also belongs to the one for partition func¬ 
tion Z\A,2^h]. The ratio between the number of such coin¬ 
cidences and the total number of simulation steps is precisely 
Z[A,2,h]/Z[hf. 

It turns out that neither of the two strategies is preferable for 
a wide range of system sizes. The thermodynamic integration 
needs a constant amount of individual simulations in order 
to densely cover the integration space, independently of the 
system size, while for the ensemble switching only a single 
simulation is needed. However, a common configuration of 
both Z\hY and Z\A, 2, h] requires the full agreement of all 
spins in part A and thus becomes a rare event when part A 
is large leading to a dramatic increase in simulation times to 
obtain small statistical errors. This drawback can be partially 
cured by a so-called increment tricli^, which subdivides the 
problem into individual simulations for subsets of a partition 
of A. 

We find that thermodynamic integration typically outper¬ 
forms the ensemble switching method for large systems (i.e. 
L > 48), while the reverse is true for smaller sizes. As such 
we employ both approaches in the following. 


Loop updates 

To efficiently sample the loop gas configurations constitut¬ 
ing the partition function of the classical toric code non-local 
update techniques are an essential requirement. While such 
non-local update techniques, e.g. in the form of loop up¬ 
dates, are well known and widely employed, for instance, in 


the context of stochastic series expansionJSl the combination 
of these techniques with the replica technique has not been 
developed so far (neither in the context of quantum systems 
nor for classical loop models). The complication one faces 
when combining these two approaches is manifesting itself at 
the boundary of the bipartition where valid loop configura¬ 
tions must be sampled at all times. In particular, this implies 
that when iteratively building up the loop update the algorithm 
needs to branch out into multiple replicas, e.g. when build¬ 
ing a loop starting from subsystem A and branching out into 
the independent subsystems Bi, B 2 , ■ ■ ■ B^, see Fig.j^for an 
illustration. One major technical achievement of this paper 
is the development of precisely such a (directed) loop update 
algorithm that efficiently samples replica-consistent loop gas 
configurations for closed loops as well as excited states with 
open loops. We refer to the appendix for a detailed descrip¬ 
tion of this update scheme and only want to briefly sketch the 
main ideas in this section. 


In an initial strategy we separate the tasks of identifying and 
accepting a new configuration. The first task is accomplished 
by performing an unbiased random walk that starts at an ar¬ 
bitrary vertex in the replicated lattice and continues through 
the lattice until it goes back to its initial vertex. If it enters 
the connected subsystem A the random walk simultaneously 
traverses in both replicas and subsequently bifurcates at the 
boundary when transitioning back into the multiple indepen¬ 
dent replicas of subsystem B. If such bifurcations occur, we 
need to not only rejoin the first vertex, but also ensure that 
all open ends at these intermediate bifurcations are eventually 
turned into closed loop configurations. Once a valid closed 
loop structure has been identified the statistical weight of the 
resulting new spin/loop configuration (in which all spins are 
flipped for all segments along the identified loop structure) is 
determined and the resulting configuration is accepted with 
Metropolis probability. This update procedure can easily be 
adapted to a case where also open loop ends are allowed (a 
scenario relevant for finite temperature simulations discussed 
in Section IV D[ ). The probability for the emergence of a new 
pair of defects is simply included in the Metropolis selection. 

A more advanced strategy to maximize the acceptance 
probability of the new configuration is to introduce statisti¬ 
cal weights already when building the update loop. To do so, 
the random walk which samples a new loop in the lattice is 
no longer unbiased but directed, i.e. the direction of turns in 
the loop update depends on the local loop configuration at ev¬ 
ery vertex and is biased with certain statistical weights. Note, 
that in our context the term directed loop refers to loops in 
the lattice and not in the imaginary-time expansion of world¬ 
line or SSE-type quantum Monte Carlo approaches. Also for 
this advanced strategy the update loop undergoes intermediate 
bifurcations if it exits subsystem A and we have to perform 
the walk until all open ends (including those occurring at the 
boundary between the two subsystems) are rejoined. Allow¬ 
ing for open loop ends (e.g. in the finite-temperature simu¬ 
lations of Section IV D| i this approach is readily modified by 
adding the possibility to stop the loop at a specific vertex and 
thereby creating an excitation. 

For further details on these two approaches, in particular 
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the detailed balance equations required to implement the di¬ 
rected loop updates, we refer to the extensive appendix of this 
manuscript. 

To obtain our simulation results we used a combination of 
both algorithms since it turns out that the unbiased random 
walk strategy is more efficient for small loop tensions h, while 
the directed walk is performing better for intermediate to large 
loop tensions h > 0.3. 

B. Mutual information 

We start our presentation of the numerical results by first 
looking at the (Renyi) mutual information - a quantity that 
tracks long range correlations between two parts of a many- 
body system. It is obtained from the Renyi entropies as 

hiA : B) = S 2 {A) + 82 ( 3 ) - S 2 {AB) . (12) 

Heuristically, the mutual information indicates how much in¬ 
formation may be gathered about subsystem A by doing all 
possible measurements in subsystem B (and vice versa). Such 
an information retrieval about part A occurs in particular in the 
presence of long-range correlations, which mediate some of 
the information about subsystem A into subsystem B where 
it can be gathered. As such, we expect the classical mutual 
information to be large in an ordered phase and easily distin¬ 
guished from the one of a disordered phase where it should 
be considerably reduced. The mutual information h as in deed 
proven to be a versatile indicator of phase transitionJ^lEl it is 
therefore natural to also consider it in the context of topologi¬ 
cal order as we will do in the following. 

Before doing do, we note that contrary to the classical en¬ 
tropy, the mutual information obeys a boundary law 

l 2 {A:B)=ce ■{£-!)+a (13) 

as all bulk contributions are canceled by definition. I 2 is in¬ 
stead sensitive to subleading contributions to the Renyi en¬ 
tropies, which we have specified here as a boundary and a 
constant contribution. Focusing on the boundary contribution 
first, note that we wrote down a term proportional to {£ — 1) 
instead of the total length i of the boundary. The reason for 
this reduced scaling form originates directly from the loop gas 
constraint that all loops need to be closed. As such only an 
even number of loop segments is allowed to cross the bound¬ 
ary between the two subsystems A and B, i.e. specifying the 
loop occupation of 1 bond segments along the boundary we 
can readily deduce the occupation of the last bond segment. 
This argument not only provides us with the £ — 1 factor, but 
also readily determines the coefficient of the boundary contri¬ 
bution to be Cf = In 2. For further details of this argument 
we refer to Ref. in which classical loop gas (and string 
net) states have been extensively discussed from a combina¬ 
torial perspective. Let us now turn to the subleading constant 
contribution a to the mutual information. This contribution 
can be readily calculated from the 0(1) contributions to the 
Renyi entropies as specified in the table of Fig. Noting 
that for a half-torus bipartition as employed in our numerical 




FIG. 6. (color online) Top panel: mutual information I 2 of the classi¬ 
cal toric code model with a magnetic field/loop tension h for varying 
system sizes. The grey shaded area indicates the intermediate regime 
he/2 < h < he- Bottom panel: Zoom of the data near he illustrating 
the crossing points for different system sizes. The inset shows a scal¬ 
ing of the crossing points with inverse system size 1/L extrapolating 
well to the expected value of he = ln(l + v/2)/2 ~ 0.440686 ... in 
the thermodynamic limit (indicated by the grey lines). 

simulations (see below) no topological 0(1) contributions are 
expected to occur, we directly probe the connectivity contri¬ 
butions to the entropies. We thus expect the constant contri¬ 
bution to the mutual information to be 

{ 0 , h < hc/n 

-;j^ln2 , hc/n < h < he 
0 , h> he, 

for arbitary Renyi entropies with n > 2. 

In our numerical setup for the classical toric code, we 
choose a symmetric bipartition of a quadratic lattice of size 
L X L with periodic boundary conditions (torus) where part 
A is precisely one half of the torus. The boundaries are two 
straight lines of length £ = L. The symmetry of this bipar¬ 
tition readily implies S 2 {A) = S 2 {B). Our results for the 
mutual information l 2 (A : B) are plotted in Fig.|^for varying 
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FIG. 7. (color online) Contributions to the finite-size scaling of the 
mutual information for the classical toric code model with magnetic 
field/loop tension h. We fit the mutual information to the form I 2 {A : 
B) = ce ■ {i — 1) + a where I is the length of the boundary of the 
bipartition. The top panel shows the boundary coefficient ci and the 
bottom panel the constant contribution a of this fit. Data is obtained 
from fitting system sizes Lmin = 12 to Lmax = 24 and Lmin = 16 
to Z/max = 24, respectively. The grey shaded area indicates the 
intermediate regime hc/2 < h < he- 


values of the loop tension h. Indeed the illustrated behav¬ 
ior of the mutual information - rescaled by the length of the 
boundary - nicely confirms our expectation; in the loop gas 
limit of small loop tension the data for different system sizes 
L collapse onto a single curve, which saturates at the expected 
value of C£ = In 2. For strong loop tension h > hf. where the 
system is in a polarized state the mutual information drops 
to zero as expected. Of interest is the intermediate regime 
/ic/2 < h < he, in which the mutual infomation first splits 
at /ic/2 for different system sizes, overshoots the value of In 2 
and subsequently exhibits a crossing points at he- The behav¬ 
ior in the intermediate region can be traced back to the fact that 
subsystem A is effectively in the polarized state, which leads 
to an overestimation of the information to be gained about this 
state from the perspective of subsystem B. 


FIG. 8. (color online) The 0(1) contribution to the Renyi entropy 
of the classical toric code model versus magnetic field/loop 
tension h determined via the Levin-Wen summation scheme. The 
0(1) contribution contains the expected universal topological con¬ 
tribution of 7 = ln(2) in the regime h < he as well as an addi¬ 
tional connectivity contribution 7 = ln(2) in the intermediate regime 
hc/2 < h < he (indicated by the grey shaded area), in which for suf¬ 
ficiently large system sizes subsystem A already transitions into the 
paramagnetic phase. 


To understand why we see a splitting versus a crossing point 
at he!2 and he, respectively, it is worthwhile taking a look at 
the subleading contributions to the mutual information. To 
this end, we have fitted our numerical results for the mutual 
information for different system sizes to the expected scaling 
behavior of Eq. 0. Results for the obtained boundary co¬ 
efficient Cf, and constant coefficient a are shown in Fig. It 
is confirmed that the boundary coefficient q attains a limit¬ 
ing curve which coincides with Fig. in the thermodynamic 
limit. More interesting is the constant contribution a whose 
sign changes explain the splitting vs. crossing points for dif¬ 
ferent system sizes. In accordance with our expectation we see 
that this constant term does follow the behavior described in 
Eq. ( [T 4 I 1 . In addition, we see that the coefficient a undergoes 
a sign change at he, which results in the crossing point of the 
mutual informatiorl^. In contrast, at he/2 it does not change 
signs but goes from zero to negative, which in turn yields the 
splitting of the mutual information (but not a crossing point). 


C. Topological entropy and connectivity contribution 


Our main quantity of interest is the classical topological en¬ 
tropy which we extract via the Levin-Wen scheme from the 
Renyi entropy. In order to reduce effects from asymmetries in 
the bipartition we choose the thickness of any building blocks 
in the scheme to be L/4. We determine this entropy for a 
range of loop tensions 0 < (1 < 0.6 since we expect it to 
reflect the phase transition at he ~ 0.44. For small system 
sizes we find a validation of the analytical expectation for the 


limits of h he and h > he provided in Sec. |IIIB| with a 
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FIG. 9. (color online) Near the critical loop tension, crossing points 
between data curves of different linear system sizes L can be iden¬ 
tified. While keeping the difference between crossing L’s constant 
(AL = 8 ) a clear tendency towards the critical point can be ob¬ 
served upon increasing L. However, an extrapolation of the location 
of the crossing points is non-trivial because a subleading logarith¬ 
mic contribution coming from corners of the subsystem is expected 
at criticality unlike in off-critical regions. 


FIG. 10. (color online) The 0(1) contribution to the Renyi entropy 
of the 2D Ising model versus inverse temperature /3 determined 
via the Levin-Wen summation scheme. While no topological con¬ 
tribution is expected, a finite 0 ( 1 ) connectivity contribution of 7 = 
ln(2) is observed in the intermediate regime /3c/2 < P < fie (indi¬ 
cated by the grey shaded area), in which for sufficiently large system 
sizes subsystem A already transitions into the low-temperature or¬ 
dered phase. 


continuous decrease of the topological entropy between the 
limiting values of In 2 and 0 around he, see Fig. However, 
upon increasing L an additional feature emerges in the range 
of /ic/2 < h < he - an interim overshooting of 7 before it 
vanishes as expected near he- This overshooting converges 
to a plateau at a value of 2 In 2 for large L. As discussed in 
Sec. IIIB| this overshooting is an artifact of the Renyi entropies 
with n > 2 and arises from the connectivity contribution. Fol¬ 
lowing the table of Fig. we see that it yields a value of 2 In 2 
for n = 2, to which the numerical data in Fig. [^nicely con¬ 
verges. 


0(1) contributions for the 2D Ising model 


The occurence of this contribution evokes the question if 
it also arises in a system without topological order. For this 
reason we consider the 2D Ising model on the square lattice 
and extract numerically the 0(1) contribution to the Renyi en¬ 
tropy when driving the Ising model through its thermal phase 


transition. Results are shown in Fig. 10 where we see the 
formation of a plateau at In 2 in the intermediate temperature 
regime between /?c/2 and j3e upon increasing system size. In 
this intermediate region, subsystem A already exhibits mag¬ 
netic order. Since here we now allow for spontaneous sym¬ 
metry breaking the two possible magnetically ordered states 
contribute a term In 2 per disconnected part of subsystem A, 
i.e. 2 In 2 for bipartition Ai and In 2 for bipartitions A 2 , A 3 , 
and A 4 . Note that this is only the case if the disconnected 
parts of subsystem A are spatially separated beyond the cor¬ 
relation length - which explains the saturation on the plateau 
only for sufficiently large system sizes. 


For higher Renyi entropies with n > 2 we thus expect a 
plateau at a level of In 2/(n — 1) between Pe/'n and Pe in the 
thermodynamic limit. 


D. Finite-temperature simulations 

For finite-temperatures, excitations that break the closed 
loop constraint of the zero-temperature loop gas are ther¬ 
mally activated. These thermally induced defects, i.e. open 
loops in the configurations, are expected to drive the system 
into a disordered phase and hence destroy its topological or¬ 
der. In fact, in the thermodynamic limit it is knowrPI that 
topological order does not survive for any finite temperature 
T > 0. For finite system sizes, however, a thermal crossover 
between the topologically ordered and the disordered phase 
takes place at some finite crossover temperature Teo{L). As 
pointed out in Ref. |35] the crossover temperature vanishes as 
Tco{L) ^ (lnL)“^, both for the quantum and classical toric 
code. 

We verify this argument by numerically calculating the van¬ 
ishing of the topological signature in the Renyi entropy via 
finite-temperature simulations, initially in the absence of a 
magnetic field, i.e. h = 0. Our results are depicted in Fig. 11 
and clearly show the expected behavior, i.e. a vanishing of 
the topological entropy for sufficiently large temperatures. To 
determine the location of the actual crossover temperature we 
identify the inflection points of the topological entropy, see 
the lower panel of Fig. [m We find that the so-determined 
crossover temperatures show a perfect 1/ In L scaling. 

At this point, one might ponder why no intermediate behav¬ 
ior appears for this transition (similar to the ones discussed 
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FIG. 11. (color online) Top panel: Finite-temperature behavior of 
topological entropy for the classical toric code model indicating the 
thermal transition from the low-temperature topologically ordered 
phase to the high-temperature paramagnet for various system sizes. 
Bottom panel: Identification of the transition temperature into the 
topologically trivial phase by determining the peak position of the 
derivative d^/dT for various system sizes. The scaling of the so- 
determined transition temperatures shows the expected logarithmic 
scaling with system size indicating that topological order is unstable 
at finite temperature in the thermodynamic hmiP. 


FIG. 12. (color online) Finite-temperature behavior of 0(1) contri¬ 
bution to the Renyi entropy of the classical toric code model in 
the presence of a finite magnetic field/loop tension h = 0.05. The 
0 ( 1 ) contribution exhibits the expected universal topological contri¬ 
bution 7 = ln( 2 ), which contributes for sufficiently large tempera¬ 
tures and vanishes at high temperatures. An additional connectivity 
contribution of 7 = ln(2) is observed in an intermediate regime. It 
can be understood from noting that the effective field is /leff = Ph 
as we introduced P in the Hamiltonian. Hence this field increases 
for r —>■ 0 and we observe the same behavior as in Fig. - read 
from right to left. The grey shaded area indicates the regime of 
he/2 < hes < he- At very low temperature even a small field 
h = 0.05 leads to the frozen trivial phase so that no topological or¬ 
der is present and thus 7 = 0 . 


can now traverse subsystem A such that the original closed- 
loop constraints for the two boundaries melt into a single con¬ 
straint for the total boundary of part A. As a consequence, not 
only the topological In 2 -contribution of bipartition ^4 van¬ 
ishes (cf. the table in Fig. but also all connectivity contri¬ 
butions (which arise from a similar boundary counting argu¬ 
ment). 


Finite magnetic field 


before) given that the replica technique again puts the system 
effectively at /3 in part B and at 2/3 in part A. A related ques¬ 
tion is whether the crossover temperature indicates whether 
subsystem A or B looses its topological order. We checked 
the latter by calculating the next higher Renyi entropy S 3 {A), 
which puts part A at temperature 3/3. As illustrated in Fig. [m 
the crossover temperatures does not move at all, thus indicat¬ 
ing that it occurs only when subsystem B undergoes order to 
the disorder transition. 

This behavior can be qualitatively understood from consid¬ 
ering, for instance, what happens to the boundary constraints 
for bipartition A 4 of the Levin-Wen summation scheme above 
the crossover temperature T^-o (where subsystem A is still 
topologically ordered while subsystem B is not). Open loops 


Finally, we turn to the combined effect of finite-temperature 
defects and loop tension in the toric code. For such a model 
the statistical weights of a loop configuration a becomes 

w{a) = exp {-pjyDa- + Phrua) , (14) 


where D^r denotes the number of loop defects and 771 ^ is 
the total magnetization. Note that in contrast to the zero- 
temperature study above, we have now introduced an explicit 
temperature also in the weight stemming from a finite loop 
tension. In particular, this implies that if the temperature is 
too low compared to the loop tension h, there is no topolog¬ 
ical order since all fluctuations of the loop gas are frozen out 
from the loop tension and the configuration is a single trivial 
state, see Fig. 12 Only if <C P~^ ^ Jy, the system sam¬ 
ples all loop configurations with approximately equal weight 
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and topological order is observed. To this extent, hes = Ph 
takes the role of an effective loop tension. Therefore, the left 
half of Fig.fT^shows precisely the same physics as Fig. [^seen 
from right to left - including the occurrence of a connectivity 
contribution for an intermediate regime hc/2 < /leff < he of 
the effective loop tension. On the other hand, for sufficiently 
high temperatures we again start to thermally activate defects 
resulting in a thermal crossover where the topological entropy 
drops to zero (see the right half of Fig. [T2) i. Again we see that 
this crossover temperature vanishes like 1/ In L with increas¬ 
ing system size - precisely mimicking the behavior seen for 
the system without magnetic field in Fig. E] 

V. DISCUSSION 

To summarize, we have studied the breakdown of topolog¬ 
ical order in classical many-body systems from a Renyi en¬ 
tropy perspective. One main result is the occurrence of an in¬ 
termediate coupling regime, in which partial subsystem order¬ 
ing leads to additional 0(1) contributions to the Renyi entropy 
- a scenario, which we extensively discussed in the context of 
the classical toric code model. 

While we have concentrated our discussion on the Levin- 
Wen summation schemd^to extract these 0(1) contributions 
it is probably worthwhile to point out that all observations 
reported here analogously apply to the alternative Kitaev- 
Preskill summation schem^ as well. We note in passing that 
in applying the Kitaev-Preskill summation scheme to detect 
topological order in classical systems, one needs to (i) adapt 
a topologically non-trivial geometry for the bipartions (e.g. 
a donut geometry) in contrast to the application for quantum 
systems and (ii) be aware of the fact that not all corner contri¬ 
butions are fully eliminated in this schemd^. 

A natural question is to ask whether such effects may also 
occur in quantum systems. To this end, it is worth pointing 
out that the replica scheme to calculate the quantum Renyi en¬ 
tropies is defined in the world-line representation of the par¬ 
tition function, i.e. replicas are stitched together in imagi¬ 
nary time or inverse temperature. Thus, the reported artifacts 
would only occur when considering thermal phase transitions 
of quantum systems and would not be expected to occur for 
zero-temperature quantum phase transitions driven by some 
coupling parameter (which is not explicitly used in the replica 
representation). 

On the other hand, we note that the classical Renyi en¬ 
tropies we discussed here are directly related^ to the basis- 
dependent Shannon entanglement entropy introduced and 
studied in Refs. [3814411 These basis-dependent Shannon and 
Renyi entropies as well as related quantities have been shown 
to exhibit unexpected additional signatures when driving a 
system through a phase transition, reflecting the same physics 
discussed here. 

Taking an even wider perspective, we might think of our 
results as indications of the limitations of the Renyi entropies 
with regard to their n = 1 counterparts, i.e. the von- 
Neumann entanglement entropy and the classical Shannon en¬ 
tropy. Other limitations of the Renyi entropies have been dis¬ 


cussed earlier, e.g. in the context of the 0( N) mod el^ as 
well as Rokhsar-Kivelson type wave functionJ^lHlEl which 
both show additional (boundary) phase transitions for Renyi 
indices n > 2. In the latter case the quantum entropy may 
be mapped onto a classical problem with different replicas 
stitched together, but all at the same temperature. Such a glu¬ 
ing can be sufficient to trigger an ordering, not in the bulk but 
at the boundary. In both cases a precise renormalization group 
analysis is needed to address the possibility of an intermedi¬ 
ate boundary phase transition. For generic, more complicated 
states such types of ordering may or may not occur - the pre¬ 
cise behavior as a function of Renyi index n would (likely) 
need to be worked out case by case depending on the model 
and on spatial dimensionality. So far all known examples refer 
to critical states without topological order in spatial dimension 
d > 1, which are driven away from criticality due to the glu¬ 
ing of several replicas. However, it appears unlikely that such 
a phenomenon could occur in gapped topologically ordered 
systems at zero temperature. Indeed in that case no connec¬ 
tivity contribution has been observed in numerical simulations 
so faP, even though these are typically limited to relatively 
small system sizes. 
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Appendix A: Replica loopgas sampling 

In this appendix we provide a detailed account of the adap¬ 
tation of non-local loop updates to the replica setup. 

The core of a Monte Carlo simulation lies in the implemen¬ 
tation of the Markov chain of configurations, i.e. the algo¬ 
rithm for selecting and accepting new configurations. This is 
sketched in the following. The major challenge when comput¬ 
ing Renyi entropies is that any configuration within the replica 
system has to satisfy the loopgas constraint in both replicas. 
Starting from a valid configuration, which could e.g. be the 
state with all spins pointing down, we will apply a modifica¬ 
tion to this state by creating or deleting loops. To do so we 
have the algorithm sample a template loop - a loop in the lat¬ 
tice which does not take care of the orientation of the spins. 
If accepted according to the detailed balance condition this 
template loop will be applied to the spins along this loop by 
flipping them. This procedure is shown in Fig. Sampling a 
template loop is particularly challenging if this loop traverses 
part A since its degrees of freedom are identified with those 
in the other replica. Thus, two additional open loop ends are 
generated in the other replica at the boundary between A and 
B. These two ends need to be linked in part B of this other 
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replica as well. Trying this, it is possible to traverse part A 
again. Obviously this procedure can easily result in a seesaw 
of linking two open loop ends in part B of either replica. 

In order to generate the template loop a first approach is to 
perform a random walk through the lattice. Once a valid loop 
is found, we apply a Metropolis probability between the cur¬ 
rent and the potential next state to decide on its acceptance. 
This is discussed in |A 1| We will then present in A 2 a more 


advanced algorithm which unifies the selection and accep¬ 
tance step in a so-called directed walk. Finally, the extension 


to finite temperatures is explained in A 3 


By a walk we mean a procedure that creates a series of adja¬ 
cent vertices whose links carry a spin. During the algorithmic 
creation of the series the walk has a head which is the the 
last appended vertex. This head “decides” which vertex (and 
thereby which spin) to append next. 



FIG. 14. (color online) Visualization of the healing process. On top, 
a template loop (green) is sampled within the first replica. It leaves 
behind two open loop ends (red circles) at the boundary between A 
and B in the other replica. The right sketch depicts the healing pro¬ 
cesses: Two random walks (blue and red) simultaneously try to meet 
each other. The dashed segments indicate parts of the template loops 
which are discarded. Internal loops (dashed red) would be allowed 
but are contracted in order to shorten the total template loop. 


1. Random walk 


The random walk is started at an arbitrary vertex of the toric 
code lattice and proceeds with equal probability of ^ to one of 
the four neighboring vertice^. Once a vertex is visited a sec¬ 
ond time, we can stop the walk and discard the first segment of 
the loop up to the first visit of this vertex. Possibly, the random 
walk has entered one or more times the connected subregion 
A so that its segments in A in the other replica also have to 
be included in the template loop. Thus, an even number of 
open template loop ends is created in the other replica at the 
boundary to part B. To match these ends, it is efficient to start 
individual random walks at every open end simultaneously. If 
the head of one of the walks hits a vertex already visited by a 
different walk, then we have created a loop that connects the 
two open ends, so that these two are healed. Typically, the 
two loops do not meet precisely at their heads so that there is 
a superfluous part between the meeting point and one of the 
heads. This part of the template loop is discarded, cf. Fig. 
14 It may also happen, that a specific random walk head hits 
a vertex twice. In this case, we can discard the resulting in¬ 
ternal loop of this walk since it does not help in linking open 
ends. Moreover, the general strategy is to keep the total length 
of the template loop as short as possible. 

In Monte Carlo parlance, the procedure described above en¬ 
sures ergodicity. What remains is to fulfill detailed balance. 
Since we have only generated a template loop so far, we can 



FIG. 13. (color online) A loop update is performed by creating a 
template loop (green) in the lattice and applying it to the spin con¬ 
figuration with its spin loops (black) if the resulting energy change is 
accepted. 


compare the total weight of the configuration prior to the ap¬ 
plication of this template loop and after it. In our special case 
we need to determine the magnetizations m-before and mafter- 
The acceptance is then decided on by the Metropolis proba¬ 
bility 

p{a —)■ cr o template loop) = min(l, (Al) 

where the function composition symbol o is used to denote the 
application of the template loop on the current configuration. 

Since this implementation of the Markov chain separates 
the selection and acceptance of new configurations we can al¬ 
most freely design the random walk procedure as long as we 
guarantee ergodicity. A shortcoming of it appears at higher 
loop tension h where the acceptance probability of an ener¬ 
getically less favourable configuration is low. In addition, for 
relatively large sizes of the subsystem A, the healing of all 
open loop ends can entail many random walks back and forth 
between the replicas. This makes the algorithm less efficient, 
especially for the parameter settings (h « 0.44, L large) we 
want to investigate to track the phase transitions. 


2. Directed walk 

A more efficient approach to the Monte Carlo update 
scheme is to unify selection and acceptance. The idea is to 
ensure detailed balance on-the-fly while generating the tem¬ 
plate loop. Strictly speaking it is no longer a template loop 
since all spins along this loop are flipped with probability 1, 
i.e. they can (and have to) be flipped directly. The walk that 
the head of the loop performs is not at random but obeys prob¬ 
ability rules depending on the value of the potential next spins 
to be visited. In general, the rule is to select the next spin 
(which unambiguously selects the next vertex) with so-called 
heat bath probability. This means, the microscopic weights Wi 
of all candidates for the next head of the loop are added up to 
a normalization constant n. In the toric code model the head 
has four possibilities to choose the next spin because bounc¬ 
ing must now be included. The only possible values for the 
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weights are Wi S {exp(h), exp(—/i), exp(2h), exp(—2/i)}, 
where the factor of two in the exponent applies in the con¬ 
nected part A which is simulated at an effective loop tension 
of 2h as discussed above. With the constant n = J2t=i 
heat bath probabilities for the four directions are defined as 


Wi W2 W3 W4 

n ^ n ^ n ' n 


(A2) 


It is not obvious to see whether this rule generates an up¬ 
date which obeys the detailed balance condition. To check it 
we need to identify a reverse update to any (completed) update 
and compare the probabilities of their occurence. If we now 
try to fractionize this update into its individual moves accord¬ 
ing to the heat bath rule and compare it to its reverse move, 
cf. Fig. 16 we see that the normalization constants thwart 


the analysis: A specific spin at position j in the loop was cho¬ 
sen to be visited (and flipped) with probability wj/rij in the 
original update. Doing its reverse, we would have the proba¬ 
bility to reflip it. The ratio of the two probabilities 

is in general not wj/wj^ as they should be from the ratio of 
the weights of the configurations that differ by flipping spin j. 
Two facts about the involved normalization constants need to 
be dealt with: (i) they originate from different configurations 
and (ii) they are shifted by a lattice position. 

Issue (i) is not a problem since the value of the normaliza¬ 
tion constant at a specific vertex at the moment when it is the 
head is the same for both the original and its reverse update. 
This can be understood from Fig. 


15 where we see that the 


coming-from and the going-to spin switch roles (and weights) 
in both situations. 






FIG. 15. (color online) Passing of the directed walk (green) at vertex 
j. The normalization constants for inverse loop update walks agree. 
In this example they both have the value of rij = 3exp(—/i) + 
exp(/i). 

The other issue (ii) can be resolved by leaving the local 
perspective from the loop head and regarding the entire loop 
update. If we consider a chain of moves and its reversion, 
we have to multiply all probabilities to obtain the total prob¬ 
ability for the walk. In the ratio of the total probabilities for 
the walk and its reverse, we realize that all normalization con¬ 
stants vanish. Since their respective numerator-denominator 
pairing is shifted by one in the product, we dub this feature of 
the microscopic rule staggered detailed balance. The stagger¬ 


ing can be seen in Fig. 16 and the ratio between a loop update 


and its reverse is given later in Eq. (A4 1 . 

We have seen that detailed balance can be satisfied in prin¬ 
ciple by the design of the walk but crucial aspects have not 


been discussed yet: How is a loop initiated and finished and 
what are the decision rules for choosing the next spin at the 
boundary between A and Bl It is clear that a loop must bifur¬ 
cate when leaving part A but how is a reunification at another 
site at the boundary accepted? The latter may fail such that 
the update is discarded for technical reasons (as opposed to 
probabilistic reasons). 

For the start of a loop update, a random spin in a random 
replica is chosen and one of its adjacent vertices is selected 
with probability of i. This spin is immediately flipped and 
the walk continues at the selected vertices following the heat 
bath rule. It is instructive to note at this point that we have 
not “paid” the flipping of this first spin in terms of acceptance 
probabilities. This will be caught up at the decision to end the 
loop. The probability to start the loop at a particular pair of a 
spin and a vertex is thus simply 


Pinit — 


1 

47V’ 
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where TV is the total number of spins in the lattice. 

In order to finish the walk it is first of all necessary that the 
loop be closed, i.e. that the head reaches precisely the other 
vertex at the initial spin — the vertex that was not chosen 
to start with. This is not sufficient, since just like anywhere 
the head is free at this vertex to chose any of the adjacent 
spins — among them the initial spin. Before we determine the 
heat bath probabilities and normalization constant we have to 
flip the initial spin once again. The loop head thus “sees” the 
original orientation of this spin before the whole loop update 
was started. Afterwards the heat bath selection is performed. 
In the case, the inital spin is chosen, it is flipped a third time 
and the loop update is successfully finished. Only now, the 
flipping of the initial spin is “paid” by a probabilistic selection 
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FIG. 16. Formation of the probabilities of two reversing update loops 
at a specific intermediate position j + 3 of the walk using heat bath 
probabilities. We compare the probabilities of selecting the encircled 
spin in the right moving loop and in its reversal left moving coun¬ 
terpart. The probability in the right moving case is Wj+a/rij+a = 
exp(/i)/(exp(—h) + 3exp(h)). Note that spins at j .. .j + 2 are 
already flipped when the loop does this selection at rij+s since we 
perform the update in-place. For the reverse move the selection prob¬ 
ability is = exp(—h)/(3exp(—/i) + exp(/i)). In par¬ 

ticular, the ratio of these to probabilities is not exp(2/i) as it should 
be if we wanted a microscopic detailed balance for the spin flip at 
i + 3. 
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according to its weight. If another but the initial is chosen, the 
loop continues and we must not forget to flip the initial spin 
again. 

We have now explained the algorithm for non-boundary¬ 


crossing loops. Fig. 17 and the following equation prove that 
the detailed balance condition is fulfilled. 
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The relative shift of the numerators again demonstrates what 
we mean by staggered detailed balance for the loop steps. The 
chronological order of the appearance of the probabilities is 
from left to right in the main numerator but from right to left 
in the main denominator. We chose to write the equation in 
this way to emphasize the reversing effect of factors that are 
written one below the other. 



FIG. 17. (color online) Completing a loop by catching up the prob¬ 
abilistic choice of the initial spin as part of the loop. The initial spin 
is drawn red in order to emphasize that it was not paid in terms of 
probabilistic weight selection until the last step. In the middle sketch 
it is temporarily flipped back since the loop head must choose it as if 
it was never flipped. 

It remains to define rules to deal with boundary crossings 
and the bifurcations of the loop that follow. In order to jus¬ 
tify the probabilities for the direction of the walk we need to 
ensure reversibility of any moves we allow. This is the rea¬ 
son why we only allow a single bifurcation of the loop and 
hope for its reunification. In case the reunification fails, the 
update is aborted and the original configuration is restored. 
We will discuss this point further at the end of this section. 
The vertices of interest for the treatment of boundary cross¬ 
ings are those which have adjacent spins in either of the bipar¬ 
titions. In the following we call them boundary vertices. In 
general we will never set a heat bath normalization constant 
at a vertex using weights from both parts, i.e. both exp(±/i) 
and exp(±2/i). In other words, a drift to or away from the 
connected part A is avoided. 

We first describe the update for the case of an initial spin 
somewhere in part B. If a boundary vertex is visited coming 
from the disconnected part B we treat all four spins as if they 
were in part B (although at least one is in part A) for the cal¬ 
culation of the heat bath weights. In other words, the head 


does not see the boundary. If the head selects a spin in part 
A, also its counterpart in the other replica is flipped. How¬ 
ever, the flipping of this counterpart is at that state not paid in 
terms of acceptance probability just as we used this wording 
above. The loop continues in part A where it uses the weights 
exp(±2/i) until it reaches again a boundary vertex. At such 
an occurrence we have to perform the mentioned bifurcation: 
First, a loop is continued in the initially chosen replica only. 
This loop should now rejoin the initial spin in the same way 
we described above for the simpler case devoid of a connected 
part A. Trying this, the loop must not enter part A again. If 
it does, we have to abort the entire update. Second, another 
loop is continued in the respective other replica. In the same 
manner, this loop is expected to rejoin the by then unpaid spin 
in part A at the entering boundary spin. To successfully finish 
the total loop update, this unification must happen without a 
prior visit of part A again. Fig. 18 illustrates the described 
procedure. 



FIG. 18. (color online) In the presence of a connected subsystem A 
loops need to bifurcate. After entering part A, all weights have to 
be taken using the effective loop tension of 2h and flipped in both 
replicas. The first spin in part A in the other at the entrance is un¬ 
paid (drawn red), just like the initial spin in part B. The bifurcation 
happens when the loop leaves part A again and must be finished in¬ 
dividually in both replicas by joining the red spins. 

A slight difference needs to be made if the initial spin of the 
update is in part A. In this case, the first visit of a boundary 
vertex leads to bifurcation and the bifurcated loop in replica 
1 of part B is free to reenter part A at any boundary vertex, 
creating there an unpaid boundary spin in replica 2. Next, the 
other bifurcated loop in replica 2 of part B must walk until it 
reaches be performed and likewise go to this unpaid bound¬ 
ary spin that its counterpart has selected for reentering part A. 
Only if this succeeds, is the subsequent loop in part A contin¬ 
ued and then the update can be completed by going back to 
the initial spin. 

Our implementation of the Monte Carlo algorithm does not 
allow multiple pairs of temporary open loops in part B in one 
of the replicas. In principle one could allow more than one 
pair of open loops and perform the loop update in the ini¬ 
tial replica until it rejoins the initial spin. Subsequently one 
could start to connect all 2no open ends in part B of the other 
replica. Doing so it would be important to start the connecting 
walks only at those Uq ends that arose from leaving part A in 
the initial replica. Special care must be taken in this approach 
to not mix the order of the start spins of the healing walks in 
part B of the other replica. They must be started in the same 
order the arose from the initial loop. This is due to the am- 
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FIG. 19. (color online) Two mutually reversing open loop updates 
including the normalization constants used for the respective choice 
of the next direction of the walk. The last normalization ua (resp. no) 
is missing and therefore will be included artificially in the decision 
to stop the loop. 


biguity of possible paths that lead to the same loop update. 
Mixing the order of starting spins would violate the detailed 
balance condition because the reversal loop of a loop could 
then be achieved through more than one move. 

We decided to allow only a single pair of bifurcations of the 
loop and avoid a long seesaw of healing processes for a single 
update step. 


3. Finite-temperature simulations 

So far, we have considered the classical tone code at zero 
temperature, i.e. the system subjected to the loopgas con¬ 
straint. By releasing this constraint and activating the vertex 
terms in the Hamilton function (j^, we can investigate the toric 
code at finite temperature. In this situation pairs of open loop 
ends are permitted at the energetical cost of two times the ver¬ 
tex coefficient J„. In practice we set Jy = 1 since only the 
ratio between Jy and h is important. A loop update in our 
Monte Carlo simulation must therefore be able to introduce 
open loops. Once again we propose two strategies. 

Our first approach consists in separating the operations of 
the loop update and the introduction of loop excitations. We 
pick two vertices of the total system at random and have the 
algorithm perform a random walk between them. At h = 0 
we only need to determine an acceptance probability for the 


open ends, which we set to Boltzmann weights exp(±J„) or 
exp(±2J„) depending on the subsystem where the open end 
is introduced. In case one of the vertices is in part A and 
the other in part B we need a third vertex in part B of the 
respective other replica and thus create a “double-tongued” 
open loop. Also at hnite loop tension we can apply Metropo¬ 
lis probabilities to the acceptance of a sampled open loop by 
balancing both the change in energy and in magnetization be¬ 
tween the current and the proposed conhguration. 


Going beyond the simple sampling technique we also 
present an algorithm that again includes the acceptance prob¬ 
ability in the selection of the subsequent configuration. We 
therefore extend our zero energy algorithm based on heat bath 
selections for the random walk at every vertex in the lattice. 
The head of the loop may take four different directions at a 
specific vertex. We add a fifth event to these four, namely 
stopping the loop and thereby creating (or annihilating) an 
open loop end in the spin conhguration. However, we will 
not include the decision of stopping the loop in the heat bath 
sum but perform it separately using the Boltzmann weight of 
an open loop end. This decision is made prior to the selection 
of the next bond to walk on. There is one issue concerning the 
detailed balance between an open loop update and its rever¬ 
sal update: Due to the staggered detailed balance philosophy 
the last heat bath normalization constant of the walk does not 
enter any of the probabilities for the loop continuation but the 
hrst one does, see Fig. 19 Roles are switched in the reversal 
update, which implies that the two normalization constants at 
the edges do not vanish in the ratio between the opposite loop 
updates. 


Our workaround is to artihcially include every heat bath 
normalization constant (and thus also the last one) into the 
decision of stopping the loop. The probability of stopping at 
the 7 ths vertex of the loop is thus Since in our model 

rij 

rij > 1 in any case, this alteration of the pure Boltzmann 
constant never leads to trivial probabilities (> 1) for stopping 
the loop. In formulas the ratio between two opposing loop 
updates is given by 
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Until now we have neglected again the bipartition and repli¬ 
cation of the lattice. Surprisingly, the situation becomes eas¬ 
ier if open loops are allowed. We can simply forbid boundary 
crossings of open loops, i.e. abort any loop update which tries 
it. The algorithm is still ergodic since an arbitrary loop up¬ 
date (including a boundary crossing) can by generated by two 
non-crossing open loops which share one of their ends at the 
boundary. Moreover, also closed loops can be created that 
way by creating an open loop whose starting and ending ver¬ 
tex coincide. This does not make the closed loop algorithm 


superfluous since it is indispensable at the hard loopgas con¬ 
straint and more efficient for very low temperature in combi¬ 
nation with the open loop algorithm. 
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